Abundant species dynamics from benthic surveys
An Bui
1/31/2022
0. set up
# source the cleaning/summarizing script, which sources the set up
source(here::here("code", "01_cleaning-and-summarizing.R"))1. Biomass
a. Guiding questions
What’s the percent community biomass that can be attributed to:
1. Laminariales, and
2. the 11 abundant species from Miller et al.?
Data question: what is the difference between the percent cover data and the percent cover that’s in the column in the biomass data?
Stray thoughts: introduced/native species trait column?
b. Test case: SCTW from 3 different sampling points
To make sure the code was working, I tested this out with “SCTW_2010-07-29”, “SCTW_2016-07-20”, and “SCTW_2020-07-29”. The issue comes from scenarios where species were counted at more than one transect - for example, at SCTW_2020-07-29, EC was at both transect 2 and 3, but I only want the dry/wet biomass for EC for the whole survey.
When I create the data frame, there are then “duplicate” observations for total species biomass for EC because it was present at both transects. To fix this, I filtered the data frame for unique observations of the total and percent biomass calculations and double checked the total dry percent to make sure everything added up to 1.
i. Calculations
This is just a test to make sure the calculations worked.
ii. Plots
These are tests of the plots.
c. Biomass data set
Since this worked, I felt comfortable doing this for the whole biomass data frame.
i. Calculations
Note: as of 2022-03-04, the code to create the prop_biomass data frame is in the 00_set-up.R script.
ii. Plots
Made some plot functions to keep things tidy.
Note: FIX THIS. YOU NEED TO REDO THESE PLOTS BECAUSE THE DATA FRAME THEY’RE BUILT ON DOESN’T TAKE THE AVERAGE AT THE SITE ANYMORE.
2. Tables of species contribution to biomass
This is a big table of percent species contribution to community biomass for each site across all sampling points. It’s useful for figuring out where to sample what. 05_species-selection-tables.R is a separate script with this code, essentially to play around with sites.
3. Urchin density at each site
Plots using stat_summary of mean urchin density at each site and a timeseries:
urchin_plot <- ggplot(urchins %>% filter(site %in% c("bull", "aque", "napl", "ivee", "mohk", "carp")), aes(x = reorder(site, density), y = density)) +
stat_summary(aes(fill = site),
fun = "mean", geom = "col") +
stat_summary(aes(group = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
labs(x = "Site", y = "Mean urchin density (number/m2)",
title = "Purple urchin density") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 35)) +
theme_bw() +
theme(legend.position = "none")
urchin_plot
urchins_time_plot <- urchins %>%
filter(site %in% c("bull", "aque", "napl", "ivee", "mohk", "carp")) %>%
filter(year > 2015) %>%
ggplot(aes(x = year, y = density)) +
# stat_summary(aes(col = site),
# fun = "sum", geom = "point", size = 3) +
# stat_summary(aes(col = site),
# fun = "sum", geom = "line") +
stat_summary(aes(col = site),
fun = "mean", geom = "point", size = 3) +
stat_summary(aes(col = site),
fun = "mean", geom = "line") +
stat_summary(aes(col = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
labs(x = "Year", y = "Sum urchin density (number/m2)",
title = "Purple urchin density through time") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(-0.1, 15)) +
theme_bw()
urchins_time_plot
Plot with pairwise comparisons of means:
# analysis of variance
urchin_aov <- aov(na.omit(urchins$density) ~ urchins$site)
# Tukey HSD
urchin_tukey <- TukeyHSD(urchin_aov)
# compact letter display
urchin_cld <- multcompLetters4(urchin_aov, urchin_tukey)
# output of compact letter display in data frame
urchin_cld_df <- enframe(urchin_cld[[1]]$Letters)
# summary data frame for plotting
urchin_plot_summary <- urchins %>%
# filter(site %in% c("bull", "aque", "ivee", "mohk", "napl", "carp")) %>%
group_by(site) %>%
summarize(mean = mean(density, na.rm = TRUE),
se = se(density)) %>%
left_join(., urchin_cld_df, by = c("site" = "name")) %>%
# mutate(value = fct_relevel(value, "a", "b", "bc", "c")) %>%
arrange(value)
# mutate(site = fct_relevel(site, "napl", "bull", "carp", "aque", "ivee", "mohk"))
# new plot
urchin_plot_new <- ggplot(urchin_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 1)) +
labs(x = "Site", y = "Mean urchin density (number/m2)",
title = "Purple urchin density") +
scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(0, 28)) +
theme_bw()
urchin_plot_new
4. Sand cover
Plot using stat_summary:
sand_plot <- ggplot(substrate %>% filter(common_name == "Sand"), aes(x = site, y = percent_cover)) +
stat_summary(aes(fill = site),
fun = "mean", geom = "col") +
stat_summary(aes(group = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
# facet_wrap(~site) +
labs(x = "Site", y = "Mean sand percent cover (%)",
title = "Sand cover across sites") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 90)) +
theme_bw() +
theme(legend.position = "none")
sand_plot
sand_time_plot <- ggplot(substrate %>% filter(common_name == "Sand"), aes(x = year, y = percent_cover)) +
stat_summary(aes(col = site),
fun = "mean", geom = "point", size = 3) +
stat_summary(aes(col = site),
fun = "mean", geom = "line") +
stat_summary(aes(col = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
# facet_wrap(~site) +
labs(x = "Site", y = "Mean sand percent cover (%)",
title = "Sand cover across sites through time") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(-1, 105)) +
theme_bw()
sand_time_plot
In case saving plots is of interest:
Plot with pairwise comparisons of means:
# filter out sand
sand <- substrate %>%
filter(common_name == "Sand")
# filter(site %in% c("bull", "aque", "ivee", "mohk", "napl", "carp"))
# analysis of variance
sand_aov <- aov(na.omit(sand$percent_cover) ~ sand$site)
# Tukey HSD
sand_tukey <- TukeyHSD(sand_aov)
# compact letter display
sand_cld <- multcompLetters4(sand_aov, sand_tukey)
# results into data frame
sand_cld_df <- enframe(sand_cld[[1]]$Letters)
# summary data frame
sand_plot_summary <- sand %>%
group_by(site) %>%
summarize(mean = mean(percent_cover, na.rm = TRUE),
se = se(percent_cover)) %>%
left_join(., sand_cld_df, by = c("site" = "name"))
# mutate(site = fct_relevel(site, "napl", "bull", "carp", "aque", "ivee", "mohk"))
# new plot
sand_plot_new <- ggplot(sand_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 1)) +
labs(x = "Site", y = "Mean sand percent cover (%)",
title = "Sand cover across sites") +
scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(0, 24)) +
theme_bw()
sand_plot_new
5. Kelp cover
Plot using stat_summary:
biomass %>%
filter(sp_code == "MAPY") %>%
filter(!(site %in% c("scdi", "sctw"))) %>%
# replace all -99999 values with 0
mutate(wm_gm2 = replace(wm_gm2, wm_gm2 < 0, 0)) %>%
ggplot(aes(x = reorder(site, wm_gm2), y = wm_gm2)) +
stat_summary(aes(fill = site),
fun = "mean", geom = "col") +
stat_summary(aes(group = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 7000)) +
theme_bw() +
theme(legend.position = "none") +
labs(x = "Site", y = "Kelp biomass (wet g/m2)",
title = "Kelp biomass across sites")
biomass %>%
filter(sp_code == "MAPY") %>%
filter(site == "carp") %>%
# filter(!(site %in% c("scdi", "sctw"))) %>%
# replace all -99999 values with 0
mutate(wm_gm2 = replace(wm_gm2, wm_gm2 < 0, 0)) %>%
ggplot(aes(x = year, y = wm_gm2)) +
stat_summary(aes(col = site),
fun = "mean", geom = "point", size = 3) +
stat_summary(aes(col = site),
fun = "mean", geom = "line") +
stat_summary(aes(col = site),
fun.data = mean_se, geom = "errorbar", width = 0.5) +
theme_bw() 
Plot with pairwise comparisons:
Note: does it make sense that mean kelp biomass is the same across all sites?
mapy <- biomass %>%
# filter(!(site %in% c("sctw", "scdi"))) %>%
filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>%
filter(sp_code == "MAPY") %>%
group_by(site, year) %>%
summarize(total = sum(dry_gm2)) ## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.mapy_aov <- aov(na.omit(mapy$total) ~ mapy$site)
mapy_tukey <- TukeyHSD(mapy_aov)
mapy_cld <- multcompLetters4(mapy_aov, mapy_tukey)
mapy_cld_df <- enframe(mapy_cld[[1]]$Letters)
mapy_plot_summary <- mapy %>%
group_by(site) %>%
summarize(mean = mean(total, na.rm = TRUE),
se = se(total)) %>%
left_join(., mapy_cld_df, by = c("site" = "name"))
mapy_plot_new <- ggplot(mapy_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 70)) +
labs(x = "Site", y = "Mean biomass (g dry/m2)",
title = "MAPY biomass across sites") +
# scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(0, 1600)) +
theme_bw()
mapy_plot_new
6. Understory biomass
macrobio <- biomass %>%
filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>%
filter(sp_code != "MAPY") %>%
group_by(site, year) %>%
summarize(total = sum(dry_gm2))## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.macrobio_aov <- aov(macrobio$total ~ macrobio$site)
macrobio_tukey <- TukeyHSD(macrobio_aov)
macrobio_cld <- multcompLetters4(macrobio_aov, macrobio_tukey)
macrobio_cld_df <- enframe(macrobio_cld[[1]]$Letters)
macrobio_plot_summary <- macrobio %>%
group_by(site) %>%
summarize(mean = mean(total, na.rm = TRUE),
se = sd(total)/sqrt(length(total))) %>%
left_join(., macrobio_cld_df, by = c("site" = "name"))
# mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))
macrobio_plot_new <- ggplot(macrobio_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 50)) +
labs(x = "Site", y = "Mean biomass (g dry/m2)",
title = "Understory macroalgal biomass across sites") +
scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(0, 1600)) +
theme_bw()
macrobio_plot_new
# all algae
all_macro <- biomass %>%
filter(site %in% c("aque", "ivee", "mohk", "napl", "carp")) %>%
group_by(site, year) %>%
summarize(total = sum(dry_gm2))## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.all_macro_aov <- aov(all_macro$total ~ all_macro$site)
all_macro_tukey <- TukeyHSD(all_macro_aov)
all_macro_cld <- multcompLetters4(all_macro_aov, all_macro_tukey)
all_macro_cld_df <- enframe(all_macro_cld[[1]]$Letters)
all_macro_plot_summary <- all_macro %>%
group_by(site) %>%
summarize(mean = mean(total, na.rm = TRUE),
se = sd(total)/sqrt(length(total))) %>%
left_join(., all_macro_cld_df, by = c("site" = "name"))
# mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))
all_macro_plot_new <- ggplot(all_macro_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 80)) +
labs(x = "Site", y = "Mean biomass (g dry/m2)",
title = "Whole community biomass across sites") +
scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(0, 4500)) +
theme_bw()
all_macro_plot_new
7. Community metrics
a. wrangling
Going to make a site by species matrix for each sampling date. There are a couple of anomalies: CARP_2000-09-08 and NAPL_2000-09-18 don’t have any observations, and AQUE_2019-07-23 is missing the sand cover data because sand cover was taken on 2019-07-24. Going to fix that later.
Note: as of 2022-03-04, the code to create community_matrix and community_metadata are in the 00_set-up.R script.
b. Species richness/diversity
i. Richness
sp_rich <- specnumber(community_matrix) %>%
enframe() %>%
left_join(., community_metadata, by = c("name" = "sample_ID")) %>%
rename(richness = value)
sp_rich_aov <- aov(sp_rich$richness ~ sp_rich$site)
sp_rich_tukey <- TukeyHSD(sp_rich_aov)
sp_rich_cld <- multcompLetters4(sp_rich_aov, sp_rich_tukey)
sp_rich_cld_df <- enframe(sp_rich_cld[[1]]$Letters)
sp_rich_plot_summary <- sp_rich %>%
group_by(site) %>%
summarize(mean = mean(richness, na.rm = TRUE),
se = sd(richness)/sqrt(length(richness))) %>%
left_join(., sp_rich_cld_df, by = c("site" = "name")) %>%
mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))
sp_rich_plot_new <- ggplot(sp_rich_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 1)) +
labs(x = "Site", y = "Mean species richness",
title = "Understory macroalgal species richness across sites") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 25)) +
theme_bw()
sp_rich_plot_new
ii. Diversity
sp_div <- diversity(community_matrix, index = "shannon") %>%
enframe() %>%
left_join(., community_metadata, by = c("name" = "sample_ID")) %>%
rename(shandiv = value)
sp_div_aov <- aov(sp_div$shandiv ~ sp_div$site)
sp_div_tukey <- TukeyHSD(sp_div_aov)
sp_div_cld <- multcompLetters4(sp_div_aov, sp_div_tukey)
sp_div_cld_df <- enframe(sp_div_cld[[1]]$Letters)
sp_div_plot_summary <- sp_div %>%
group_by(site) %>%
summarize(mean = mean(shandiv, na.rm = TRUE),
se = sd(shandiv)/sqrt(length(shandiv))) %>%
left_join(., sp_div_cld_df, by = c("site" = "name")) %>%
mutate(site = fct_relevel(site, "napl", "carp", "aque", "ivee", "mohk"))
sp_div_plot_new <- ggplot(sp_div_plot_summary, aes(x = site, y = mean)) +
geom_col() +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.5) +
geom_text(aes(label = value, y = mean + se + 0.2)) +
labs(x = "Site", y = "Mean Shannon diversity",
title = "Understory macroalgal Shannon diversity across sites") +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 2.75)) +
theme_bw()
sp_div_plot_new
sp_div %>%
ggplot(aes(x = site, y = shandiv)) +
geom_violin() +
geom_point() +
geom_text(data = sp_div_cld_df, aes(x = name, y = 3, label = value))
c. SIMPER analysis
Wanted to know which species significantly contributed to differences in community composition.
sim <- with(community_metadata, simper(community_matrix, site))
head(summary(sim))## $carp_napl
## average sd ratio ava avb cumsum
## PTCA 2.455694e-01 0.2556539579 0.96055407 8.694819937 8.264586e+01 0.2861201
## CC 8.088812e-02 0.1036097668 0.78069975 5.306398810 2.057173e+01 0.3803653
## DL 7.458567e-02 0.1265995199 0.58914656 10.447767857 1.363259e+01 0.4672672
## CYOS 4.770110e-02 0.0789746402 0.60400527 8.626758988 5.942733e+00 0.5228451
## CO 4.735428e-02 0.0928932663 0.50977087 14.833928571 6.656250e-01 0.5780190
## R 4.551648e-02 0.0492601324 0.92400233 2.605208333 1.165536e+01 0.6310516
## EC 4.407481e-02 0.0601931871 0.73222256 11.176785714 6.742411e+00 0.6824044
## POLA 3.749029e-02 0.0957137185 0.39169195 8.316964286 2.410714e-02 0.7260855
## GR 3.347978e-02 0.0866076882 0.38656827 10.279017857 8.537946e-01 0.7650937
## RAT 3.026258e-02 0.0548837476 0.55139412 1.575297619 4.980134e+00 0.8003536
## LAFA 2.548957e-02 0.0530884206 0.48013421 0.010628472 6.305970e+00 0.8300522
## LS 1.959313e-02 0.0410853960 0.47688803 2.070833333 1.560417e+00 0.8528807
## DP 1.532958e-02 0.0345931828 0.44313886 4.162500000 4.901786e-01 0.8707417
## TALE 1.341341e-02 0.0489963948 0.27376316 2.595535714 1.607143e-02 0.8863700
## FB 1.161466e-02 0.0245496910 0.47310826 1.654910714 4.738839e-01 0.8999026
## CP 8.888653e-03 0.0327503077 0.27140671 0.843750000 7.955357e-01 0.9102590
## BR 7.159473e-03 0.0127488668 0.56157720 0.295238095 1.559226e+00 0.9186007
## Nandersoniana 6.707379e-03 0.0227786354 0.29445920 0.894940476 1.414583e+00 0.9264157
## BO 6.499400e-03 0.0204183675 0.31831142 0.901785714 8.416667e-01 0.9339883
## SCCA 6.424783e-03 0.0191308503 0.33583365 0.552455357 6.026786e-01 0.9414741
## EGME 5.748187e-03 0.0398427455 0.14427185 0.027107639 1.855239e+00 0.9481714
## CAL 5.699457e-03 0.0128044920 0.44511390 0.120238095 1.127232e+00 0.9548120
## CRYP 4.006483e-03 0.0314669032 0.12732373 2.049702381 0.000000e+00 0.9594801
## BF 3.767757e-03 0.0207867151 0.18125794 0.635119048 2.886905e-01 0.9638700
## DU 3.638079e-03 0.0166719754 0.21821525 1.478571429 0.000000e+00 0.9681089
## BLD 3.435456e-03 0.0093860603 0.36601687 0.483333333 3.475397e-01 0.9721116
## AMZO 3.360020e-03 0.0144763468 0.23210417 1.112202381 0.000000e+00 0.9760265
## COF 3.099095e-03 0.0105681556 0.29324841 0.000000000 6.026786e-01 0.9796373
## ER 2.247823e-03 0.0056742276 0.39614611 0.387946429 1.448661e-01 0.9822563
## SAMU 1.902505e-03 0.0255191202 0.07455214 0.007921627 3.156746e-01 0.9844730
## CG 1.500490e-03 0.0048250891 0.31097671 0.029910714 3.290179e-01 0.9862213
## CF 1.421864e-03 0.0031479900 0.45167369 0.072916667 3.541667e-01 0.9878779
## UEC 1.301632e-03 0.0128091952 0.10161704 0.269419643 1.584821e-02 0.9893945
## GS 1.118705e-03 0.0055895448 0.20014235 0.100446429 1.674107e-01 0.9906979
## GYSP 8.693371e-04 0.0032248178 0.26957711 0.054464286 1.543155e-01 0.9917108
## UV 6.940577e-04 0.0030300964 0.22905465 0.144642857 1.607143e-02 0.9925195
## HAGL 6.369483e-04 0.0051652848 0.12331330 0.000000000 1.446429e-01 0.9932616
## PL 6.322902e-04 0.0037412891 0.16900330 0.138318452 6.287202e-02 0.9939983
## FR 5.986664e-04 0.0033116541 0.18077565 0.107142857 1.071429e-01 0.9946958
## PRSP 5.896792e-04 0.0029515596 0.19978563 0.088020833 8.802083e-02 0.9953829
## AU 5.667135e-04 0.0038155426 0.14852763 0.171875000 0.000000e+00 0.9960432
## BRA 5.341714e-04 0.0023441884 0.22787050 0.063541667 8.169643e-02 0.9966655
## STIN 5.236085e-04 0.0030633004 0.17092954 0.118005952 7.261905e-02 0.9972756
## CZ 4.396347e-04 0.0029117167 0.15098816 0.138318452 0.000000e+00 0.9977878
## UBB 3.847832e-04 0.0028032770 0.13726194 0.064583333 0.000000e+00 0.9982362
## LI 3.363351e-04 0.0028926203 0.11627351 0.110937500 1.584821e-02 0.9986280
## LX 2.976988e-04 0.0029800372 0.09989769 0.043750000 1.458333e-02 0.9989749
## BPSE 2.271476e-04 0.0021710678 0.10462482 0.000000000 4.821429e-02 0.9992396
## SAHO 1.373665e-04 0.0007833223 0.17536392 0.035117064 8.591268e-04 0.9993996
## FG 1.170322e-04 0.0010029952 0.11668269 0.011904762 1.190476e-02 0.9995360
## FTHR 1.082206e-04 0.0008847429 0.12231873 0.000000000 2.083333e-02 0.9996621
## SAFU 9.881098e-05 0.0008503173 0.11620483 0.031250000 5.208333e-03 0.9997772
## EA 9.535253e-05 0.0009729069 0.09800787 0.000000000 1.468304e-02 0.9998883
## IR 7.439040e-05 0.0010242725 0.07262754 0.000000000 1.257440e-02 0.9999750
## PHSE 2.149273e-05 0.0002862812 0.07507559 0.009077381 0.000000e+00 1.0000000
## FASP 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## NEO 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## SELO 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
##
## $carp_bull
## average sd ratio ava avb cumsum
## PTCA 3.090448e-01 0.2289500398 1.34983527 8.694819937 135.07343809 0.3373353
## BF 9.197659e-02 0.0935287221 0.98340477 0.635119048 38.49206349 0.4377316
## CYOS 7.641745e-02 0.0858260143 0.89037628 8.626758988 24.61623939 0.5211444
## DL 5.730248e-02 0.0861525975 0.66512771 10.447767857 16.36428571 0.5836924
## CO 4.502185e-02 0.0635774250 0.70814207 14.833928571 8.28333333 0.6328356
## Nandersoniana 3.806181e-02 0.0390967152 0.97352958 0.894940476 15.31984127 0.6743817
## EC 3.673142e-02 0.0396696416 0.92593280 11.176785714 0.51071429 0.7144756
## EGME 3.428148e-02 0.0471466249 0.72712479 0.027107639 12.71089601 0.7518952
## CC 3.417775e-02 0.0385564501 0.88643409 5.306398810 10.56250000 0.7892017
## POLA 2.826186e-02 0.0632361257 0.44692592 8.316964286 2.25000000 0.8200507
## GS 2.662457e-02 0.0483545079 0.55061199 0.100446429 9.01785714 0.8491125
## GR 2.208290e-02 0.0626752735 0.35233834 10.279017857 0.26785714 0.8732169
## R 1.505678e-02 0.0173127522 0.86969305 2.605208333 4.35714286 0.8896520
## AU 9.361261e-03 0.0259640187 0.36054747 0.171875000 1.77777778 0.8998702
## DP 9.322010e-03 0.0237906499 0.39183502 4.162500000 0.00000000 0.9100456
## TALE 8.699694e-03 0.0303837829 0.28632689 2.595535714 0.19285714 0.9195416
## FB 6.980982e-03 0.0144868999 0.48188240 1.654910714 0.43869048 0.9271617
## BO 6.928598e-03 0.0125652895 0.55140776 0.901785714 2.24444444 0.9347245
## LS 6.914296e-03 0.0170950689 0.40446141 2.070833333 0.00000000 0.9422718
## CRYP 6.521753e-03 0.0269508422 0.24198697 2.049702381 1.38571429 0.9493905
## RAT 5.641507e-03 0.0075841288 0.74385696 1.575297619 1.10337302 0.9555485
## BR 4.059314e-03 0.0077867227 0.52131231 0.295238095 1.40238095 0.9599794
## CF 3.273523e-03 0.0053257990 0.61465387 0.072916667 1.44444444 0.9635526
## BLD 3.057056e-03 0.0073184387 0.41771977 0.483333333 0.77026455 0.9668895
## DU 2.733843e-03 0.0129779054 0.21065366 1.478571429 0.00000000 0.9698736
## UV 2.558582e-03 0.0068386079 0.37413782 0.144642857 0.57857143 0.9726664
## PHSE 2.372076e-03 0.0052753723 0.44965083 0.009077381 0.87142857 0.9752556
## AMZO 2.338950e-03 0.0104913894 0.22294000 1.112202381 0.00000000 0.9778087
## SCCA 2.247923e-03 0.0065271360 0.34439647 0.552455357 0.17857143 0.9802624
## PL 2.243585e-03 0.0047384041 0.47348967 0.138318452 0.63710317 0.9827113
## CP 2.092718e-03 0.0080118151 0.26120399 0.843750000 0.00000000 0.9849956
## FTHR 1.915062e-03 0.0035805577 0.53485030 0.000000000 0.75000000 0.9870860
## NEO 1.836746e-03 0.0064514187 0.28470421 0.000000000 0.62500000 0.9890909
## CAL 1.222971e-03 0.0049922565 0.24497362 0.120238095 0.24047619 0.9904258
## ER 1.160133e-03 0.0031800399 0.36481718 0.387946429 0.02619048 0.9916921
## BRA 1.114865e-03 0.0027286734 0.40857386 0.063541667 0.38730159 0.9929090
## STIN 9.685879e-04 0.0028830848 0.33595540 0.118005952 0.31468254 0.9939663
## UEC 8.021044e-04 0.0081596284 0.09830158 0.269419643 0.04226190 0.9948418
## LX 7.337854e-04 0.0028315173 0.25914919 0.043750000 0.23333333 0.9956428
## PRSP 6.526429e-04 0.0029732107 0.21950779 0.088020833 0.20119048 0.9963552
## LAFA 6.249676e-04 0.0021643085 0.28876088 0.010628472 0.16287632 0.9970373
## CZ 5.641034e-04 0.0026609633 0.21199219 0.138318452 0.06706349 0.9976531
## SAFU 5.601568e-04 0.0016520556 0.33906653 0.031250000 0.13888889 0.9982645
## GYSP 4.853155e-04 0.0013990578 0.34688738 0.054464286 0.14523810 0.9987943
## FR 3.065458e-04 0.0019088689 0.16059029 0.107142857 0.04761905 0.9991289
## LI 2.271339e-04 0.0021488350 0.10570097 0.110937500 0.00000000 0.9993768
## UBB 2.179603e-04 0.0015475531 0.14084192 0.064583333 0.00000000 0.9996147
## IR 1.151722e-04 0.0009342512 0.12327751 0.000000000 0.03353175 0.9997404
## CG 8.870658e-05 0.0004796712 0.18493204 0.029910714 0.01329365 0.9998373
## SAHO 8.669998e-05 0.0005131600 0.16895310 0.035117064 0.00000000 0.9999319
## FG 3.518608e-05 0.0004979940 0.07065563 0.011904762 0.00000000 0.9999703
## SAMU 2.721045e-05 0.0003949935 0.06888834 0.007921627 0.00000000 1.0000000
## BPSE 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
## COF 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
## EA 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
## FASP 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
## HAGL 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
## SELO 0.000000e+00 0.0000000000 NaN 0.000000000 0.00000000 1.0000000
##
## $carp_mohk
## average sd ratio ava avb cumsum
## PTCA 1.629987e-01 0.2123082228 0.76774558 8.694819937 44.0066146 0.1798691
## CO 9.507108e-02 0.1112406399 0.85464345 14.833928571 13.4456250 0.2847801
## CYOS 9.466860e-02 0.1104015005 0.85749376 8.626758988 13.9713325 0.3892470
## EC 8.792094e-02 0.1173889353 0.74897125 11.176785714 0.4331250 0.4862678
## R 7.908317e-02 0.0763753306 1.03545443 2.605208333 17.1562500 0.5735361
## CC 6.670531e-02 0.0640267406 1.04183512 5.306398810 13.7840625 0.6471454
## GR 4.800132e-02 0.1003171248 0.47849582 10.279017857 2.3906250 0.7001149
## POLA 4.388456e-02 0.1093719281 0.40124151 8.316964286 0.1012500 0.7485415
## DL 4.053225e-02 0.1070211051 0.37873137 10.447767857 0.8156250 0.7932689
## RAT 2.072370e-02 0.0241639399 0.85762927 1.575297619 4.1456250 0.8161375
## DP 1.640438e-02 0.0368369271 0.44532423 4.162500000 0.4387500 0.8342397
## TALE 1.632516e-02 0.0564391536 0.28925237 2.595535714 0.2025000 0.8522546
## LS 1.526779e-02 0.0401103140 0.38064501 2.070833333 0.0000000 0.8691026
## FB 1.430644e-02 0.0337190026 0.42428425 1.654910714 0.3196875 0.8848897
## BLD 1.271000e-02 0.0277770116 0.45757261 0.483333333 2.2765000 0.8989152
## BO 1.139346e-02 0.0266472330 0.42756637 0.901785714 1.5150000 0.9114879
## EGME 1.078418e-02 0.0243393811 0.44307554 0.027107639 2.1278385 0.9233883
## BR 7.688069e-03 0.0178998054 0.42950575 0.295238095 1.3175000 0.9318721
## GS 5.697721e-03 0.0151997199 0.37485696 0.100446429 1.2656250 0.9381595
## CRYP 4.895620e-03 0.0334951620 0.14615903 2.049702381 0.1212500 0.9435618
## DU 4.324915e-03 0.0179240565 0.24129109 1.478571429 0.0675000 0.9483344
## SCCA 4.130474e-03 0.0133412827 0.30960093 0.552455357 0.1406250 0.9528924
## GYSP 3.949196e-03 0.0101914712 0.38750005 0.054464286 0.6100000 0.9572503
## AMZO 3.806647e-03 0.0161119745 0.23626201 1.112202381 0.0000000 0.9614509
## Nandersoniana 3.702192e-03 0.0153584594 0.24105232 0.894940476 0.4850000 0.9655363
## CP 3.662915e-03 0.0142207179 0.25757595 0.843750000 0.0000000 0.9695783
## AU 3.542542e-03 0.0120498855 0.29398967 0.171875000 0.5468750 0.9734875
## CZ 3.366317e-03 0.0101240971 0.33250538 0.138318452 0.4753125 0.9772023
## PRSP 2.803459e-03 0.0090408161 0.31008917 0.088020833 0.3696875 0.9802959
## ER 2.278642e-03 0.0069553852 0.32760833 0.387946429 0.0000000 0.9828104
## BF 1.917702e-03 0.0123580476 0.15517841 0.635119048 0.0000000 0.9849265
## BRA 1.724404e-03 0.0051040139 0.33785261 0.063541667 0.3050000 0.9868294
## LAFA 1.590329e-03 0.0053810939 0.29554016 0.010628472 0.2301229 0.9885844
## PL 1.539362e-03 0.0068096210 0.22605690 0.138318452 0.2112500 0.9902830
## CAL 1.452365e-03 0.0056392407 0.25754623 0.120238095 0.1262500 0.9918857
## UEC 1.311896e-03 0.0142724374 0.09191814 0.269419643 0.0000000 0.9933334
## FR 1.215001e-03 0.0047748041 0.25446099 0.107142857 0.3250000 0.9946742
## LI 1.040688e-03 0.0048743876 0.21350122 0.110937500 0.1331250 0.9958226
## UV 9.769952e-04 0.0037999680 0.25710618 0.144642857 0.0337500 0.9969007
## UBB 7.340540e-04 0.0035975763 0.20404126 0.064583333 0.0775000 0.9977107
## STIN 4.335038e-04 0.0024839272 0.17452356 0.118005952 0.0381250 0.9981891
## CF 4.089273e-04 0.0012163657 0.33618776 0.072916667 0.0656250 0.9986403
## PHSE 3.729136e-04 0.0024715066 0.15088513 0.009077381 0.0381250 0.9990518
## SAFU 2.418915e-04 0.0013145763 0.18400720 0.031250000 0.0437500 0.9993188
## CG 2.110499e-04 0.0009967199 0.21174442 0.029910714 0.0209375 0.9995517
## SAHO 1.532170e-04 0.0008932640 0.17152489 0.035117064 0.0000000 0.9997207
## LX 1.291713e-04 0.0012608312 0.10244928 0.043750000 0.0000000 0.9998633
## FG 6.468797e-05 0.0009041190 0.07154807 0.011904762 0.0000000 0.9999347
## SAMU 5.921697e-05 0.0008700064 0.06806499 0.007921627 0.0000000 1.0000000
## BPSE 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## COF 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## EA 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## FASP 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## FTHR 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## HAGL 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## IR 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## NEO 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
## SELO 0.000000e+00 0.0000000000 NaN 0.000000000 0.0000000 1.0000000
##
## $carp_ivee
## average sd ratio ava avb cumsum
## DL 0.1312661303 0.1900446250 0.69071214 10.447767857 2.624500e+01 0.1470734
## PTCA 0.1197036199 0.1761482632 0.67956174 8.694819937 1.998016e+01 0.2811918
## CYOS 0.0865615525 0.1105684061 0.78287782 8.626758988 1.390637e+01 0.3781772
## CO 0.0820668440 0.1242198862 0.66065786 14.833928571 7.484583e+00 0.4701267
## EC 0.0752750791 0.0890583815 0.84523296 11.176785714 1.705000e+00 0.5544665
## POLA 0.0615202171 0.1206187433 0.51003862 8.316964286 4.230000e+00 0.6233950
## GR 0.0335516730 0.0962484522 0.34859442 10.279017857 3.125000e-02 0.6609870
## BF 0.0323666054 0.0709527129 0.45617150 0.635119048 5.496667e+00 0.6972513
## GS 0.0299702168 0.0615556399 0.48688011 0.100446429 6.625000e+00 0.7308305
## CC 0.0286887442 0.0535456148 0.53578139 5.306398810 3.966806e+00 0.7629740
## R 0.0237801211 0.0444181865 0.53536902 2.605208333 3.117778e+00 0.7896178
## Nandersoniana 0.0225936897 0.0484140311 0.46667648 0.894940476 4.149444e+00 0.8149322
## TALE 0.0180237829 0.0547605608 0.32913803 2.595535714 5.550000e-01 0.8351265
## FB 0.0177747700 0.0339181192 0.52404940 1.654910714 1.228333e+00 0.8550417
## DP 0.0158521388 0.0366420429 0.43262159 4.162500000 1.650000e-01 0.8728028
## LS 0.0154228949 0.0362885040 0.42500774 2.070833333 1.088889e-01 0.8900829
## RAT 0.0108499994 0.0170140440 0.63770844 1.575297619 7.351389e-01 0.9022395
## EGME 0.0099991100 0.0339286238 0.29471016 0.027107639 2.511308e+00 0.9134427
## BR 0.0076581706 0.0204210909 0.37501281 0.295238095 9.300000e-01 0.9220231
## BO 0.0076310186 0.0240526797 0.31726272 0.901785714 7.294444e-01 0.9305730
## SAMU 0.0061392105 0.0494033702 0.12426704 0.007921627 2.982759e+00 0.9374515
## AU 0.0058244785 0.0133882138 0.43504523 0.171875000 1.234722e+00 0.9439774
## CRYP 0.0055297062 0.0352598296 0.15682737 2.049702381 1.077778e-01 0.9501730
## PHSE 0.0048949395 0.0184521368 0.26527765 0.009077381 7.794444e-01 0.9556574
## SCCA 0.0048623369 0.0135937524 0.35768909 0.552455357 2.812500e-01 0.9611052
## BLD 0.0046276268 0.0126573174 0.36560881 0.483333333 3.888148e-01 0.9662901
## DU 0.0040411007 0.0180548791 0.22382319 1.478571429 0.000000e+00 0.9708179
## AMZO 0.0039123107 0.0158838362 0.24630767 1.112202381 2.805556e-02 0.9752013
## CP 0.0036636905 0.0137773343 0.26592158 0.843750000 0.000000e+00 0.9793062
## CF 0.0026185647 0.0057238008 0.45748704 0.072916667 4.277778e-01 0.9822401
## ER 0.0024485328 0.0063627871 0.38482079 0.387946429 5.958333e-02 0.9849835
## BRA 0.0018300085 0.0083293562 0.21970588 0.063541667 2.033333e-01 0.9870338
## UV 0.0015325063 0.0059160217 0.25904339 0.144642857 1.500000e-01 0.9887509
## FTHR 0.0015199490 0.0051650240 0.29427724 0.000000000 3.013889e-01 0.9904539
## UEC 0.0014030451 0.0140551331 0.09982439 0.269419643 2.958333e-02 0.9920259
## UBB 0.0013207516 0.0060098354 0.21976503 0.064583333 2.238889e-01 0.9935057
## GYSP 0.0010254566 0.0039650690 0.25862264 0.054464286 1.355556e-01 0.9946546
## CAL 0.0006641022 0.0041075407 0.16167879 0.120238095 0.000000e+00 0.9953987
## FR 0.0006278653 0.0033488718 0.18748563 0.107142857 6.666667e-02 0.9961022
## CZ 0.0004967955 0.0032102137 0.15475465 0.138318452 0.000000e+00 0.9966588
## SAFU 0.0004838790 0.0020137357 0.24028924 0.031250000 8.750000e-02 0.9972009
## LX 0.0004374759 0.0025369251 0.17244338 0.043750000 5.444444e-02 0.9976911
## LI 0.0004015853 0.0031774971 0.12638417 0.110937500 2.958333e-02 0.9981410
## PL 0.0003666459 0.0023285160 0.15745905 0.138318452 0.000000e+00 0.9985518
## PRSP 0.0003335074 0.0026474804 0.12597162 0.088020833 2.347222e-02 0.9989255
## LAFA 0.0003151463 0.0019141828 0.16463750 0.010628472 4.222917e-02 0.9992786
## STIN 0.0003001699 0.0023514742 0.12765179 0.118005952 0.000000e+00 0.9996149
## SAHO 0.0001533313 0.0008698014 0.17628314 0.035117064 2.222222e-04 0.9997867
## CG 0.0001254836 0.0007117883 0.17629349 0.029910714 9.305556e-03 0.9999273
## FG 0.0000648734 0.0008859051 0.07322839 0.011904762 0.000000e+00 1.0000000
## BPSE 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## COF 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## EA 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## FASP 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## HAGL 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## IR 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## NEO 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## SELO 0.0000000000 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
##
## $carp_aque
## average sd ratio ava avb cumsum
## PTCA 1.505860e-01 0.1974155792 0.76278691 8.694819937 3.010669e+01 0.1697838
## DL 1.487785e-01 0.2109229092 0.70536925 10.447767857 3.133125e+01 0.3375296
## EC 9.373259e-02 0.1318370641 0.71097302 11.176785714 2.414547e+00 0.4432119
## CYOS 7.683385e-02 0.1115759669 0.68862363 8.626758988 9.861216e+00 0.5298411
## CO 6.821917e-02 0.1113665372 0.61256436 14.833928571 3.213362e+00 0.6067573
## POLA 6.812189e-02 0.1443652138 0.47187189 8.316964286 4.922845e+00 0.6835638
## GR 3.566193e-02 0.1042796397 0.34198366 10.279017857 0.000000e+00 0.7237722
## R 3.239649e-02 0.0460650028 0.70327774 2.605208333 4.272629e+00 0.7602988
## CC 2.486131e-02 0.0525042771 0.47351023 5.306398810 2.185345e+00 0.7883296
## TALE 2.096720e-02 0.0646640638 0.32424806 2.595535714 7.913793e-01 0.8119699
## LS 1.855623e-02 0.0477778732 0.38838533 2.070833333 6.336207e-02 0.8328918
## DP 1.797162e-02 0.0411647354 0.43657796 4.162500000 1.512931e-01 0.8531545
## FB 1.778685e-02 0.0395719638 0.44948109 1.654910714 5.298491e-01 0.8732090
## RAT 1.310815e-02 0.0230978745 0.56750477 1.575297619 8.663793e-01 0.8879882
## AU 1.104456e-02 0.0285461994 0.38690112 0.171875000 1.689655e+00 0.9004408
## GS 9.445987e-03 0.0247688436 0.38136569 0.100446429 1.575970e+00 0.9110910
## BO 9.042056e-03 0.0267112618 0.33851100 0.901785714 7.400862e-01 0.9212858
## BR 6.343935e-03 0.0141916302 0.44701944 0.295238095 9.219828e-01 0.9284386
## Nandersoniana 5.999693e-03 0.0200053150 0.29990497 0.894940476 7.943966e-01 0.9352031
## CRYP 5.872691e-03 0.0362177102 0.16214971 2.049702381 1.254310e-01 0.9418245
## CP 5.463515e-03 0.0189751600 0.28792985 0.843750000 2.443966e-01 0.9479846
## SCCA 4.829829e-03 0.0153858522 0.31391362 0.552455357 1.454741e-01 0.9534301
## BF 4.710111e-03 0.0159120540 0.29600897 0.635119048 5.853448e-01 0.9587407
## DU 4.279185e-03 0.0189715645 0.22555786 1.478571429 0.000000e+00 0.9635654
## AMZO 4.217860e-03 0.0176383907 0.23912950 1.112202381 0.000000e+00 0.9683210
## BLD 3.594995e-03 0.0112036398 0.32087741 0.483333333 2.250000e-01 0.9723743
## ER 2.716476e-03 0.0082045793 0.33109263 0.387946429 2.489224e-02 0.9754371
## UV 2.512290e-03 0.0068189010 0.36843037 0.144642857 3.840517e-01 0.9782697
## EGME 1.845034e-03 0.0128303836 0.14380196 0.027107639 2.172195e-01 0.9803499
## BRA 1.681693e-03 0.0072122150 0.23317283 0.063541667 2.103448e-01 0.9822460
## UEC 1.657701e-03 0.0160432854 0.10332680 0.269419643 2.295259e-02 0.9841151
## SAMU 1.499332e-03 0.0128068709 0.11707249 0.007921627 1.950359e-01 0.9858055
## CF 1.467704e-03 0.0038851556 0.37777215 0.072916667 2.338362e-01 0.9874604
## GYSP 1.200532e-03 0.0049878660 0.24069055 0.054464286 1.709052e-01 0.9888139
## SELO 1.167799e-03 0.0090619374 0.12886856 0.000000000 1.396552e-01 0.9901306
## CAL 1.043564e-03 0.0054220778 0.19246564 0.120238095 4.353448e-02 0.9913072
## PHSE 9.191940e-04 0.0037901100 0.24252437 0.009077381 2.366379e-01 0.9923436
## UBB 8.769062e-04 0.0045678448 0.19197373 0.064583333 8.017241e-02 0.9933323
## STIN 7.196498e-04 0.0055606242 0.12941888 0.118005952 3.943966e-02 0.9941437
## LI 5.673895e-04 0.0040027860 0.14174865 0.110937500 4.590517e-02 0.9947834
## FASP 5.374982e-04 0.0033386566 0.16099236 0.000000000 1.051724e-01 0.9953894
## CZ 5.341173e-04 0.0035014174 0.15254316 0.138318452 0.000000e+00 0.9959917
## FR 5.130059e-04 0.0031828609 0.16117760 0.107142857 5.172414e-02 0.9965701
## LX 4.781302e-04 0.0026632471 0.17952904 0.043750000 6.336207e-02 0.9971092
## NEO 4.775126e-04 0.0031700373 0.15063312 0.000000000 1.454741e-01 0.9976475
## PRSP 4.326203e-04 0.0029773153 0.14530552 0.088020833 3.642241e-02 0.9981353
## PL 4.299539e-04 0.0024933571 0.17243975 0.138318452 1.821121e-02 0.9986201
## FTHR 3.774483e-04 0.0016700850 0.22600542 0.000000000 8.297414e-02 0.9990456
## LAFA 1.859816e-04 0.0015872062 0.11717544 0.010628472 1.303376e-02 0.9992553
## SAHO 1.732297e-04 0.0010044809 0.17245696 0.035117064 8.620690e-05 0.9994507
## IR 1.413228e-04 0.0016593587 0.08516712 0.000000000 1.821121e-02 0.9996100
## CG 1.088585e-04 0.0007335932 0.14839084 0.029910714 0.000000e+00 0.9997327
## FG 1.084103e-04 0.0010890267 0.09954784 0.011904762 8.620690e-03 0.9998550
## SAFU 9.846824e-05 0.0009588661 0.10269237 0.031250000 0.000000e+00 0.9999660
## EA 3.016904e-05 0.0003365065 0.08965367 0.000000000 7.212644e-03 1.0000000
## BPSE 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## COF 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
## HAGL 0.000000e+00 0.0000000000 NaN 0.000000000 0.000000e+00 1.0000000
##
## $napl_bull
## average sd ratio ava avb cumsum
## PTCA 2.776579e-01 2.080389e-01 1.33464424 8.264586e+01 135.07343809 0.3617611
## BF 8.074493e-02 8.432750e-02 0.95751603 2.886905e-01 38.49206349 0.4669639
## CYOS 5.971607e-02 7.076090e-02 0.84391334 5.942733e+00 24.61623939 0.5447680
## DL 5.356836e-02 7.445622e-02 0.71946116 1.363259e+01 16.36428571 0.6145624
## CC 4.729296e-02 5.755267e-02 0.82173349 2.057173e+01 10.56250000 0.6761805
## Nandersoniana 3.318287e-02 3.419358e-02 0.97044149 1.414583e+00 15.31984127 0.7194145
## EGME 3.210429e-02 4.576318e-02 0.70153107 1.855239e+00 12.71089601 0.7612433
## R 2.422633e-02 2.789667e-02 0.86843085 1.165536e+01 4.35714286 0.7928078
## GS 2.277136e-02 4.222268e-02 0.53931592 1.674107e-01 9.01785714 0.8224767
## CO 2.027309e-02 2.699765e-02 0.75092060 6.656250e-01 8.28333333 0.8488905
## EC 1.705371e-02 1.890369e-02 0.90213624 6.742411e+00 0.51071429 0.8711098
## LAFA 1.436343e-02 2.842505e-02 0.50530894 6.305970e+00 0.16287632 0.8898240
## RAT 1.319429e-02 2.137444e-02 0.61729274 4.980134e+00 1.10337302 0.9070149
## AU 7.111229e-03 1.983495e-02 0.35852019 0.000000e+00 1.77777778 0.9162801
## POLA 6.086456e-03 2.467037e-02 0.24671121 2.410714e-02 2.25000000 0.9242102
## BR 5.410701e-03 7.936182e-03 0.68177624 1.559226e+00 1.40238095 0.9312598
## BO 5.390379e-03 7.358635e-03 0.73252434 8.416667e-01 2.24444444 0.9382829
## LS 4.775916e-03 1.503448e-02 0.31766418 1.560417e+00 0.00000000 0.9445055
## CAL 3.269553e-03 7.043646e-03 0.46418477 1.127232e+00 0.24047619 0.9487654
## CF 3.112419e-03 4.764849e-03 0.65320408 3.541667e-01 1.44444444 0.9528205
## CRYP 3.107866e-03 8.153521e-03 0.38116855 0.000000e+00 1.38571429 0.9568698
## FB 3.038833e-03 9.027028e-03 0.33663711 4.738839e-01 0.43869048 0.9608291
## CP 2.635013e-03 1.416004e-02 0.18608804 7.955357e-01 0.00000000 0.9642623
## GR 2.506314e-03 6.053939e-03 0.41399716 8.537946e-01 0.26785714 0.9675277
## BLD 2.356245e-03 5.192892e-03 0.45374423 3.475397e-01 0.77026455 0.9705977
## SCCA 2.227789e-03 7.582121e-03 0.29382132 6.026786e-01 0.17857143 0.9735003
## PHSE 2.040969e-03 4.585201e-03 0.44512093 0.000000e+00 0.87142857 0.9761595
## UV 1.867275e-03 5.411747e-03 0.34504100 1.607143e-02 0.57857143 0.9785924
## PL 1.853173e-03 4.108275e-03 0.45108307 6.287202e-02 0.63710317 0.9810069
## FTHR 1.686795e-03 3.162811e-03 0.53332164 2.083333e-02 0.75000000 0.9832046
## NEO 1.565331e-03 5.568965e-03 0.28108120 0.000000e+00 0.62500000 0.9852441
## COF 1.533881e-03 5.049379e-03 0.30377609 6.026786e-01 0.00000000 0.9872426
## DP 1.393048e-03 5.969179e-03 0.23337351 4.901786e-01 0.00000000 0.9890576
## BRA 9.744933e-04 2.345415e-03 0.41548855 8.169643e-02 0.38730159 0.9903272
## SAMU 9.207575e-04 1.299543e-02 0.07085243 3.156746e-01 0.00000000 0.9915269
## TALE 8.511697e-04 4.662908e-03 0.18254054 1.607143e-02 0.19285714 0.9926359
## STIN 8.238719e-04 2.514940e-03 0.32759101 7.261905e-02 0.31468254 0.9937093
## CG 8.067520e-04 2.563852e-03 0.31466398 3.290179e-01 0.01329365 0.9947604
## GYSP 6.895111e-04 1.830906e-03 0.37659561 1.543155e-01 0.14523810 0.9956588
## LX 6.181498e-04 2.548389e-03 0.24256496 1.458333e-02 0.23333333 0.9964642
## PRSP 6.033651e-04 2.373910e-03 0.25416514 8.802083e-02 0.20119048 0.9972503
## SAFU 4.267935e-04 1.321708e-03 0.32291070 5.208333e-03 0.13888889 0.9978064
## ER 3.595011e-04 9.671849e-04 0.37169846 1.448661e-01 0.02619048 0.9982748
## HAGL 3.488862e-04 2.849984e-03 0.12241691 1.446429e-01 0.00000000 0.9987293
## FR 2.792231e-04 1.556732e-03 0.17936493 1.071429e-01 0.04761905 0.9990931
## CZ 2.122696e-04 1.437574e-03 0.14765819 0.000000e+00 0.06706349 0.9993697
## IR 1.332323e-04 9.488075e-04 0.14042080 1.257440e-02 0.03353175 0.9995433
## UEC 1.267877e-04 1.072501e-03 0.11821686 1.584821e-02 0.04226190 0.9997085
## BPSE 1.245917e-04 1.214452e-03 0.10259088 4.821429e-02 0.00000000 0.9998708
## EA 4.426948e-05 4.468248e-04 0.09907569 1.468304e-02 0.00000000 0.9999285
## FG 3.185882e-05 3.203910e-04 0.09943731 1.190476e-02 0.00000000 0.9999700
## LI 2.051270e-05 2.690408e-04 0.07624384 1.584821e-02 0.00000000 0.9999967
## SAHO 2.515191e-06 3.552632e-05 0.07079796 8.591268e-04 0.00000000 1.0000000
## AMZO 0.000000e+00 0.000000e+00 NaN 0.000000e+00 0.00000000 1.0000000
## DU 0.000000e+00 0.000000e+00 NaN 0.000000e+00 0.00000000 1.0000000
## FASP 0.000000e+00 0.000000e+00 NaN 0.000000e+00 0.00000000 1.0000000
## SELO 0.000000e+00 0.000000e+00 NaN 0.000000e+00 0.00000000 1.0000000
## UBB 0.000000e+00 0.000000e+00 NaN 0.000000e+00 0.00000000 1.0000000d. NMDS
i. Jaccard (presence/absence)
community_mds_jacc <- metaMDS(community_presabs, distance = "jaccard")## Run 0 stress 0.2324413
## Run 1 stress 0.2333624
## Run 2 stress 0.2500419
## Run 3 stress 0.2528057
## Run 4 stress 0.2347918
## Run 5 stress 0.2339213
## Run 6 stress 0.2327236
## ... Procrustes: rmse 0.006726982 max resid 0.09394191
## Run 7 stress 0.2350828
## Run 8 stress 0.2372256
## Run 9 stress 0.2339764
## Run 10 stress 0.2334229
## Run 11 stress 0.2368917
## Run 12 stress 0.2389732
## Run 13 stress 0.2364392
## Run 14 stress 0.2327612
## ... Procrustes: rmse 0.006176838 max resid 0.07131762
## Run 15 stress 0.2418717
## Run 16 stress 0.240909
## Run 17 stress 0.2332244
## Run 18 stress 0.2495368
## Run 19 stress 0.2332678
## Run 20 stress 0.2381245
## *** No convergence -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 13: stress ratio > sratmax
## 3: scale factor of the gradient < sfgrmincommunity_dist_jacc <- vegdist(community_presabs, method = "jaccard")
stressplot(community_mds_jacc)
community_mds_sites_jacc <- scores(community_mds_jacc, display = "sites") %>%
as.data.frame() %>%
rownames_to_column("sample_ID") %>%
left_join(., community_metadata, by = "sample_ID")
community_mds_spp_jacc <- scores(community_mds_jacc, display = "species") %>%
as.data.frame() %>%
rownames_to_column("sp_code") %>%
left_join(., coarse_traits, by = "sp_code") %>%
filter(sp_code %in% algae_proposal)
perm <- adonis2(community_dist_jacc ~ site, data = community_metadata)
perm## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = community_dist_jacc ~ site, data = community_metadata)
## Df SumOfSqs R2 F Pr(>F)
## site 5 29.442 0.15199 22.906 0.001 ***
## Residual 639 164.265 0.84801
## Total 644 193.707 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1disper <- betadisper(community_dist_jacc, community_metadata$site)
anova(disper)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 5 1.2762 0.255236 32.14 < 2.2e-16 ***
## Residuals 639 5.0745 0.007941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(disper)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## bull-aque -0.106271846 -0.146136993 -0.066406699 0.0000000
## carp-aque 0.015459985 -0.015289770 0.046209740 0.7042879
## ivee-aque 0.012470074 -0.023310680 0.048250828 0.9191448
## mohk-aque -0.053008215 -0.099713866 -0.006302564 0.0156226
## napl-aque -0.071037639 -0.101787394 -0.040287884 0.0000000
## carp-bull 0.121731831 0.084100706 0.159362956 0.0000000
## ivee-bull 0.118741920 0.076899161 0.160584679 0.0000000
## mohk-bull 0.053263631 0.001766330 0.104760932 0.0378176
## napl-bull 0.035234207 -0.002396918 0.072865332 0.0814327
## ivee-carp -0.002989911 -0.036263533 0.030283711 0.9998482
## mohk-carp -0.068468200 -0.113282143 -0.023654257 0.0002122
## napl-carp -0.086497624 -0.114290051 -0.058705196 0.0000000
## mohk-ivee -0.065478289 -0.113882879 -0.017073698 0.0016896
## napl-ivee -0.083507713 -0.116781335 -0.050234091 0.0000000
## napl-mohk -0.018029424 -0.062843367 0.026784519 0.8601005fit_jacc <- community_metadata %>%
select(sample_ID, site, urchin_density, mean_sand, mean_daily_umol, kelp_dry) %>%
mutate(across(.cols = 3:4, ~replace(., is.na(.), 0))) %>%
envfit(community_mds_jacc, ., perm = 999, na.rm = TRUE)
fit_vect_jacc <- scores(fit, display = "vectors") %>%
as.data.frame()
nmds_plot_jacc <- ggplot() +
# x and y axes
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
# points for sites
geom_point(data = community_mds_sites_jacc, aes(x = NMDS1, y = NMDS2, color = site)) +
stat_ellipse(data = community_mds_sites_jacc, aes(x = NMDS1, y = NMDS2, color = site), level = 0.95, size = 1) +
# # arrows for species
# geom_segment(data = community_mds_spp_jacc, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm")), color = "blue") +
# geom_text(data = community_mds_spp_jacc, aes(x = NMDS1, y = NMDS2, label = scientific_name), color = "blue") +
# # arrows for environmental characteristics (vector fitting)
geom_segment(data = fit_vect_jacc, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
geom_text(data = fit_vect_jacc, aes(x = NMDS1, y = NMDS2, label = rownames(fit_vect_jacc))) +
# theme
theme_bw() +
guides(colour = guide_legend(nrow = 5)) +
theme(legend.position = "bottom", legend.box = "vertical", legend.margin = margin())
nmds_plot_jacc
ii. Bray-Curtis
community_mds_bray <- metaMDS(community_matrix, distance = "bray")## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2397652
## Run 1 stress 0.2460827
## Run 2 stress 0.2590958
## Run 3 stress 0.243352
## Run 4 stress 0.2495171
## Run 5 stress 0.2458121
## Run 6 stress 0.2492311
## Run 7 stress 0.2471866
## Run 8 stress 0.2478991
## Run 9 stress 0.24517
## Run 10 stress 0.2592116
## Run 11 stress 0.2403553
## Run 12 stress 0.2508534
## Run 13 stress 0.2502277
## Run 14 stress 0.2489769
## Run 15 stress 0.249233
## Run 16 stress 0.2485628
## Run 17 stress 0.2489904
## Run 18 stress 0.2502947
## Run 19 stress 0.2506252
## Run 20 stress 0.247575
## *** No convergence -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmaxcommunity_dist_bray <- vegdist(community_matrix, method = "bray")
stressplot(community_mds_bray)
community_mds_sites_bray <- scores(community_mds_bray, display = "sites") %>%
as.data.frame() %>%
rownames_to_column("sample_ID") %>%
left_join(., community_metadata, by = "sample_ID")
community_mds_spp_bray <- scores(community_mds_bray, display = "species") %>%
as.data.frame() %>%
rownames_to_column("sp_code") %>%
left_join(., coarse_traits, by = "sp_code") %>%
filter(sp_code %in% algae_proposal)
perm <- adonis2(community_dist_bray ~ site, data = community_metadata)
perm## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = community_dist_bray ~ site, data = community_metadata)
## Df SumOfSqs R2 F Pr(>F)
## site 5 36.368 0.15729 23.853 0.001 ***
## Residual 639 194.856 0.84271
## Total 644 231.225 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1disper <- betadisper(community_dist_bray, community_metadata$site)
anova(disper)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 5 1.4579 0.291571 25.107 < 2.2e-16 ***
## Residuals 639 7.4208 0.011613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(disper)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## bull-aque -0.166971431 -0.21517981 -0.118763050 0.0000000
## carp-aque -0.028129109 -0.06531437 0.009056153 0.2570883
## ivee-aque -0.007283007 -0.05055219 0.035986173 0.9968071
## mohk-aque -0.077425080 -0.13390559 -0.020944570 0.0013792
## napl-aque -0.065291455 -0.10247672 -0.028106194 0.0000100
## carp-bull 0.138842322 0.09333551 0.184349131 0.0000000
## ivee-bull 0.159688425 0.10908854 0.210288305 0.0000000
## mohk-bull 0.089546351 0.02727136 0.151821338 0.0006349
## napl-bull 0.101679976 0.05617317 0.147186784 0.0000000
## ivee-carp 0.020846102 -0.01939124 0.061083442 0.6766880
## mohk-carp -0.049295972 -0.10348886 0.004896921 0.0986720
## napl-carp -0.037162347 -0.07077135 -0.003553341 0.0203913
## mohk-ivee -0.070142074 -0.12867709 -0.011607059 0.0085374
## napl-ivee -0.058008449 -0.09824579 -0.017771109 0.0006074
## napl-mohk 0.012133625 -0.04205927 0.066326518 0.9879368fit <- community_metadata %>%
select(sample_ID, site, urchin_density, mean_sand, mean_daily_umol, kelp_dry) %>%
mutate(across(.cols = 3:4, ~replace(., is.na(.), 0))) %>%
envfit(community_mds_bray, ., perm = 999, na.rm = TRUE)
fit_vect_bray <- scores(fit, display = "vectors") %>%
as.data.frame()
nmds_plot_bray <- ggplot() +
geom_point(data = community_mds_sites_bray, aes(x = NMDS1, y = NMDS2, color = site, shape = site)) +
# scale_color_manual(values = pal) +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
stat_ellipse(data = community_mds_sites_bray, aes(x = NMDS1, y = NMDS2, color = site), level = 0.95, size = 1) +
geom_segment(data = community_mds_spp_bray, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
geom_text(data = community_mds_spp_bray, aes(x = NMDS1, y = NMDS2, label = scientific_name)) +
# geom_segment(data = fit_vect_bray, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.2, "cm"))) +
# geom_text(data = fit_vect_bray, aes(x = NMDS1, y = NMDS2, label = rownames(fit_vect_bray))) +
theme_bw() +
guides(colour = guide_legend(nrow = 5)) +
theme(legend.position = "bottom", legend.box = "vertical", legend.margin = margin())
nmds_plot_bray
8. Trait calculations
a. proportion of traits
“Mean proportion of biomass” for each species per trait across sampling dates: matrix where row is sampling date, and columns are traits spread across, and cells are total biomass that can be “assigned” to that trait. The prop_biomass data frame is already set up.
i. growth form
Going to test this out with a couple of sites just to make sure the math is right.
Now that it seems ok, going to do the rest of the data set.
prop_traits <- prop_biomass %>%
filter(site %in% c("aque", "napl", "ivee", "mohk", "carp")) %>%
select(sample_ID, site, year, sp_code, growth_form,
dry_gm2, wm_gm2, total_dry, total_wet, percent_sp_dry, percent_sp_wet) %>%
# create a new column where total biomass for the whole survey is calculated
# group_by(site, year) %>%
# mutate(total_dry = sum(dry_gm2),
# total_wet = sum(wm_gm2)) %>%
# ungroup() %>%
# calculate biomass within each growth form
group_by(sample_ID, growth_form) %>%
mutate(total_gf_dry = sum(dry_gm2),
total_gf_wet = sum(wm_gm2)) %>%
ungroup() %>%
select(sample_ID, growth_form, total_gf_dry, total_dry, total_gf_wet, total_wet) %>%
unique() %>%
group_by(sample_ID, growth_form) %>%
summarize(percent_gf_dry = total_gf_dry/total_dry,
percent_gf_wet = total_gf_wet/total_wet) %>%
ungroup() %>%
# double check: should sum to 1
# group_by(sample_ID) %>%
# mutate(test = sum(unique(percent_gf_dry)))
# make the data frame wider
select(sample_ID, growth_form, percent_gf_dry) %>%
unique() %>%
pivot_wider(names_from = "growth_form", values_from = "percent_gf_dry")## `summarise()` has grouped output by 'sample_ID'. You can override using the `.groups` argument.GLM???
9. Community-weighted means
These are calculated using the functcomp function in the FD package. The function requires a matrix with functional traits and a matrix of species abundances (or presence-absence). In the trait matrix, rows are species and columns are traits. In the species matrix, rows are samples and columns are species.
The first step is to get the coarse_traits data frame into a functional form (lol).
trait_matrix <- coarse_traits %>%
select(sp_code, growth_morph, growth_form, pigment_type,
height_class, longevity, posture, branching_yn) %>%
filter(sp_code != "MAPY") %>%
column_to_rownames("sp_code")
# making sure the rownames in the trait matrix == columns in the community matrix
rownames(trait_matrix) == colnames(community_matrix)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [22] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEThen the community matrix has to be a matrix, not a data frame (why?).
comm_matrix <- as.matrix(community_matrix)Then, you’re ready to calculate community-weighted means.
CWM_cat <- functcomp(trait_matrix, comm_matrix,
CWM.type = c("all")) %>%
rownames_to_column("sample_ID") %>%
left_join(., community_metadata, by = "sample_ID") %>%
# to get the glmm to work!!!!!!
mutate(across(.cols = everything(), .fns = ~ replace(., . == 1, 0.999)))
ggplot(CWM_cat, aes(x = year, y = mean_daily_umol)) +
geom_point(aes(col = site)) +
geom_smooth(method = "lm", aes(col = site))## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 262 rows containing non-finite values (stat_smooth).## Warning: Removed 262 rows containing missing values (geom_point).
For categorical traits, the output is the exact same as the prop_traits data frame up top - each cell is the proportion of the community biomass that is attributed to a specific category (for example, within growth_morph, at CARP in 2000, there is 92% solitary and 8 % aggregated).
Now to plot?
a. GLMs
Kill me!
lm_plot <- ggplot(CWM_cat,
aes(x = mean_sand, y = growth_form_leathery_macrophyte)) +
# scale_x_continuous(expand = c(0, 0)) +
# scale_y_continuous(expand = c(0, 0)) +
# coord_cartesian(ylim = c(-0.001, 1)) +
geom_point(aes(col = site)) +
geom_smooth(method = "lm", se = FALSE, aes(col = site)) +
labs(x = "Mean sand cover", y = "Proportion biomass (dry g/m2)",
title = "Proportion leathery macrophytes")
# missing values for mean sand
lm_plot## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
ggplot(CWM_cat, aes(x = growth_form_leathery_macrophyte)) +
geom_histogram(binwidth = 0.1)
model_lm_01a <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluationmodel_lm_01b <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluationmodel_lm_01c <- glmmTMB(growth_form_leathery_macrophyte ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
data = CWM_cat, ziformula = ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry, family = beta_family(link = "logit"))## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN function evaluation# plot residuals
model_lm_01c_resid <- simulateResiduals(model_lm_01c)
plot(model_lm_01c_resid)
model_lm_nozi_resid <- simulateResiduals(model_lm_nozi)
plot(model_lm_nozi_resid)
# look at the model
model_lm_01c## Formula: growth_form_leathery_macrophyte ~ urchin_density + mean_sand +
## mean_daily_umol + kelp_dry + site + (1 | year)
## Zero inflation: ~urchin_density + mean_sand + mean_daily_umol + kelp_dry
## Data: CWM_cat
## AIC BIC logLik df.resid
## 50.04870 113.09144 -9.02435 364
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## year (Intercept) 0.2156
##
## Number of obs: 380 / Conditional model: year, 14
##
## Dispersion parameter for beta family (): 3.17
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) urchin_density mean_sand mean_daily_umol kelp_dry sitecarp
## -1.096e+00 -9.588e-04 1.950e-02 -2.229e-06 5.067e-04 -4.798e-01
## siteivee sitemohk sitenapl
## -1.938e-01 1.267e+00 1.425e+00
##
## Zero-inflation model:
## (Intercept) urchin_density mean_sand mean_daily_umol kelp_dry
## -2.309e+00 9.233e-02 -1.312e-01 -4.245e-05 -1.509e-03summary(model_lm_01c)## Family: beta ( logit )
## Formula: growth_form_leathery_macrophyte ~ urchin_density + mean_sand +
## mean_daily_umol + kelp_dry + site + (1 | year)
## Zero inflation: ~urchin_density + mean_sand + mean_daily_umol + kelp_dry
## Data: CWM_cat
##
## AIC BIC logLik deviance df.resid
## 50.0 113.1 -9.0 18.0 364
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.04646 0.2156
## Number of obs: 380, groups: year, 14
##
## Dispersion parameter for beta family (): 3.17
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.096e+00 2.119e-01 -5.175 2.28e-07 ***
## urchin_density -9.588e-04 7.810e-03 -0.123 0.902296
## mean_sand 1.950e-02 4.554e-03 4.282 1.85e-05 ***
## mean_daily_umol -2.229e-06 1.298e-05 -0.172 0.863629
## kelp_dry 5.067e-04 1.506e-04 3.365 0.000766 ***
## sitecarp -4.798e-01 2.173e-01 -2.208 0.027230 *
## siteivee -1.938e-01 2.142e-01 -0.905 0.365540
## sitemohk 1.267e+00 2.405e-01 5.270 1.37e-07 ***
## sitenapl 1.425e+00 2.143e-01 6.650 2.92e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.309e+00 3.801e-01 -6.073 1.25e-09 ***
## urchin_density 9.233e-02 1.577e-02 5.855 4.76e-09 ***
## mean_sand -1.312e-01 9.638e-02 -1.361 0.174
## mean_daily_umol -4.245e-05 3.670e-05 -1.157 0.247
## kelp_dry -1.509e-03 7.035e-04 -2.145 0.032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# significant effects of sand and kelp biomass, not of urchin density and year
# why is it an increasing slope? this doesn't make sense
AIC(model_lm_01a, model_lm_01b, model_lm_01c)## df AIC
## model_lm_01a 13 66.42352
## model_lm_01b 12 188.27592
## model_lm_01c 16 50.04870# plot predicted values on top of predictor
temp <- ggpredict(model_lm, terms = "urchin_density")
# predict(model_lm)
plot(temp, add.data = TRUE)
b. Checking correlations
Checking to see if urchins, sand, kelp, and irradiance are at all correlated.
ggplot(CWM_cat, aes(x = urchin_density, y = mean_sand)) +
geom_point(aes(col = site)) +
geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
# create correlation matrix
cor_mat <- cor(CWM_cat[, c("urchin_density", "mean_sand", "mean_daily_umol", "kelp_dry")], use = "na.or.complete")
# plot correlation matrix
corrplot(cor_mat, type = "lower", addCoef.col = "black")
# 1) sand and urchins are slightly negatively correlated,
# 2) irradiance is moderately positively correlated with urchin density,
# 3) urchin density is moderately negatively correlated with kelp biomass,
# 4) sand cover is slightly negatively correlated with kelp biomass, and
# 5) average daily irradiance is slightly negatively correlated with kelp
# check for collinearity in the model
check_collinearity(model_lm_01c)## # Check for Multicollinearity
##
## * conditional component:
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## urchin_density 1.34 1.16 0.74
## mean_sand 1.53 1.24 0.65
## mean_daily_umol 1.06 1.03 0.95
## kelp_dry 1.62 1.27 0.62
## site 2.27 1.51 0.44
##
## * zero inflated component:
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## urchin_density 1.16 1.08 0.86
## mean_sand 1.02 1.01 0.98
## mean_daily_umol 1.03 1.01 0.97
## kelp_dry 1.12 1.06 0.89# plot VIF for each variable (if >5, inflates)
plot(check_collinearity(model_lm_01c))
10. CWM for categorical and continuous traits
a. fixed CWM using FD::dbFD()
Use community_matrix and tbspp_matrix from set up script.
cwm_commat <- community_matrix %>%
select(c(rownames(tbspp_matrix))) %>%
# take out 0 sites
mutate(total = rowSums(across(where(is.numeric)))) %>%
filter(total > 0) %>%
select(-total) %>%
as.matrix()
# making sure the rownames in the trait matrix == columns in the community matrix
rownames(tbspp_matrix) == colnames(cwm_commat)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUECWM_fixed <- functcomp(x = tbspp_matrix, a = cwm_commat,
CWM.type = c("all")) %>%
rownames_to_column("sample_ID") %>%
rename_with(.fn = ~ paste0(., "_fixed", sep = ""), .cols = 2:28) %>%
left_join(., community_metadata, by = "sample_ID")i. Thickness
# plot
ggplot(CWM_fixed, aes(x = urchin_density, y = thickness_mm_fixed)) +
geom_point(aes(col = site)) +
geom_smooth(aes(col = site), method = "lm") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
ggplot(CWM_fixed, aes(x = thickness_mm_fixed)) +
geom_histogram(binwidth = 0.05)
model_thickness_01a <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
model_thickness_01b <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
data = CWM_fixed, family = gaussian(link = "identity")) ## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small eigenvalues detected. See
## vignette('troubleshooting')model_thickness_01c <- glmmTMB(thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
AIC(model_thickness_01a, model_thickness_01b, model_thickness_01c)## df AIC
## model_thickness_01a 11 -803.6142
## model_thickness_01b 7 -733.3296
## model_thickness_01c 8 -791.2190# plot residuals
model_thickness_01a_resid <- simulateResiduals(model_thickness_01a)
plot(model_thickness_01a_resid)
model_thickness_01b_resid <- simulateResiduals(model_thickness_01b)
plot(model_thickness_01b_resid)
model_thickness_02_resid <- simulateResiduals(model_thickness_02)
plot(model_thickness_02_resid)
summary(model_thickness_01)## Family: gaussian ( identity )
## Formula: thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + site + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## -800.2 -758.4 411.1 -822.2 321
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.000482 0.02195
## Residual 0.004674 0.06837
## Number of obs: 332, groups: year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 0.00467
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.136e-01 1.571e-02 19.966 < 2e-16 ***
## urchin_density -1.847e-03 5.515e-04 -3.349 0.000811 ***
## mean_sand 5.136e-04 3.068e-04 1.674 0.094111 .
## mean_daily_umol -1.751e-06 9.921e-07 -1.765 0.077553 .
## kelp_dry 2.230e-05 1.050e-05 2.123 0.033746 *
## sitecarp 8.434e-02 1.485e-02 5.679 1.36e-08 ***
## siteivee 7.621e-03 1.457e-02 0.523 0.601022
## sitemohk 6.433e-02 1.654e-02 3.889 0.000100 ***
## sitenapl 5.450e-02 1.459e-02 3.736 0.000187 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_thickness_02)## Linear mixed model fit by REML ['lmerMod']
## Formula: thickness_mm_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
##
## REML criterion at convergence: -724.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1667 -0.7238 -0.0136 0.6721 3.3539
##
## Random effects:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.0005743 0.02396
## site (Intercept) 0.0012382 0.03519
## Residual 0.0047831 0.06916
## Number of obs: 332, groups: year, 14; site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.546e-01 1.939e-02 18.292
## urchin_density -1.727e-03 5.435e-04 -3.177
## mean_sand 4.464e-04 3.065e-04 1.457
## mean_daily_umol -1.770e-06 9.560e-07 -1.851
## kelp_dry 2.299e-05 1.061e-05 2.166
##
## Correlation of Fixed Effects:
## (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.217
## mean_sand -0.171 0.075
## mean_dly_ml -0.235 -0.041 0.103
## kelp_dry -0.250 0.122 0.135 -0.146
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# check collinearity
check_collinearity(model_thickness_01)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## urchin_density 1.20 1.10 0.83
## mean_sand 1.56 1.25 0.64
## mean_daily_umol 1.06 1.03 0.94
## kelp_dry 1.64 1.28 0.61
## site 2.36 1.54 0.42plot(check_collinearity(model_thickness_01))
ii. LDMC/TDMC
ggplot(CWM_fixed, aes(x = urchin_density, y = tdmc_fixed)) +
geom_point(aes(col = site)) +
geom_smooth(aes(col = site), method = "lm")## `geom_smooth()` using formula 'y ~ x'
ggplot(CWM_fixed, aes(x = tdmc_fixed)) +
geom_histogram(binwidth = 0.01)
model_tdmc_01 <- glmmTMB(tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
testDispersion(model_tdmc_01)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sidedtestResiduals(model_tdmc_01)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.208, p-value = 1.335e-12
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate expectations
##
## data: simulationOutput
## outliers at both margin(s) = 2, observations = 324, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0007484357 0.0221193534
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.00617284## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.208, p-value = 1.335e-12
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.1352, p-value = 0.512
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate expectations
##
## data: simulationOutput
## outliers at both margin(s) = 2, observations = 324, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0007484357 0.0221193534
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.00617284model_tdmc_02 <- lmer(tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed)## Warning: Some predictor variables are on very different scales: consider rescaling# plot residuals
model_tdmc_01_resid <- simulateResiduals(model_tdmc_01)
plot(model_tdmc_01_resid)## qu = 0.75, log(sigma) = -2.428833 : outer Newton did not converge fully.
model_tdmc_02_resid <- simulateResiduals(model_tdmc_02)
plot(model_tdmc_02_resid)## qu = 0.75, log(sigma) = -1.775149 : outer Newton did not converge fully.
summary(model_tdmc_01)## Family: gaussian ( identity )
## Formula: tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +
## (1 | site) + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## -674.0 -643.8 345.0 -690.0 316
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.003701 0.06084
## year (Intercept) 0.000248 0.01575
## Residual 0.006423 0.08015
## Number of obs: 324, groups: site, 5; year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 0.00642
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.498e-01 3.013e-02 8.292 < 2e-16 ***
## urchin_density -1.221e-03 6.529e-04 -1.870 0.0615 .
## mean_sand -6.485e-04 3.582e-04 -1.811 0.0702 .
## mean_daily_umol -5.158e-06 1.196e-06 -4.313 1.61e-05 ***
## kelp_dry -1.119e-05 1.231e-05 -0.909 0.3631
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_tdmc_02)## Linear mixed model fit by REML ['lmerMod']
## Formula: tdmc_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
##
## REML criterion at convergence: -611.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85301 -0.55790 -0.05674 0.41836 2.68137
##
## Random effects:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.000293 0.01712
## site (Intercept) 0.004645 0.06816
## Residual 0.006486 0.08053
## Number of obs: 324, groups: year, 14; site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.493e-01 3.261e-02 7.645
## urchin_density -1.203e-03 6.296e-04 -1.910
## mean_sand -6.431e-04 3.581e-04 -1.796
## mean_daily_umol -5.095e-06 1.036e-06 -4.919
## kelp_dry -1.138e-05 1.240e-05 -0.918
##
## Correlation of Fixed Effects:
## (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.133
## mean_sand -0.117 0.051
## mean_dly_ml -0.148 -0.122 0.131
## kelp_dry -0.184 0.136 0.125 -0.113
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingiii. SA:P ratio
ggplot(CWM_fixed, aes(x = urchin_density, y = ratio_fixed)) +
geom_point(aes(col = site)) +
geom_smooth(aes(col = site), method = "lm")## `geom_smooth()` using formula 'y ~ x'
ggplot(CWM_fixed, aes(x = ratio_fixed)) +
geom_histogram(binwidth = 1)
# full
model_ratio_01 <- glmmTMB(ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
# plot residuals
model_ratio_01_resid <- simulateResiduals(model_ratio_01)
plot(model_ratio_01_resid)
model_ratio_02_resid <- simulateResiduals(model_ratio_02)
plot(model_ratio_02_resid)
summary(model_ratio_01)## Family: gaussian ( identity )
## Formula: ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## 2016.8 2047.1 -1000.4 2000.8 316
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 16.084 4.011
## year (Intercept) 3.826 1.956
## Residual 25.062 5.006
## Number of obs: 324, groups: site, 5; year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 25.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.209e+00 2.004e+00 3.597 0.000321 ***
## urchin_density 1.377e-01 4.130e-02 3.333 0.000858 ***
## mean_sand 9.882e-02 2.245e-02 4.401 1.08e-05 ***
## mean_daily_umol 3.901e-04 7.413e-05 5.262 1.42e-07 ***
## kelp_dry 8.832e-04 7.924e-04 1.115 0.265049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_ratio_02)## Linear mixed model fit by REML ['lmerMod']
## Formula: ratio_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
##
## REML criterion at convergence: 2061.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2541 -0.5628 -0.1867 0.5733 3.2969
##
## Random effects:
## Groups Name Variance Std.Dev.
## year (Intercept) 1.915 1.384
## site (Intercept) 28.577 5.346
## Residual 25.241 5.024
## Number of obs: 329, groups: year, 14; site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.9180720 2.5116439 2.754
## urchin_density 0.1225312 0.0397446 3.083
## mean_sand 0.1017156 0.0224186 4.537
## mean_daily_umol 0.0004240 0.0000678 6.253
## kelp_dry 0.0009980 0.0007681 1.299
##
## Correlation of Fixed Effects:
## (Intr) urchn_ mn_snd mn_dl_
## urchn_dnsty -0.112
## mean_sand -0.096 0.053
## mean_dly_ml -0.129 -0.086 0.124
## kelp_dry -0.146 0.122 0.131 -0.121
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingAIC(model_ratio_01, model_ratio_02)## Warning in AIC.default(model_ratio_01, model_ratio_02): models are not all fitted to the same number of
## observations## df AIC
## model_ratio_01 8 2016.829
## model_ratio_02 8 2077.731iv. FvFm
ggplot(CWM_fixed, aes(x = urchin_density, y = fvfm_fixed)) +
geom_point(aes(col = site)) +
geom_smooth(aes(col = site), method = "lm")## `geom_smooth()` using formula 'y ~ x'
ggplot(CWM_fixed, aes(x = fvfm_fixed)) +
geom_histogram(binwidth = 0.01)
# full
model_fvfm_01a <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed, family = Gamma(link = "logit"))
# site as fixed instead of random effect
model_fvfm_01b <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
data = CWM_fixed, family = Gamma(link = "logit"))
# no site
model_fvfm_01c <- glmmTMB(fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
data = CWM_fixed, family = Gamma(link = "logit"))
# plot residuals
model_fvfm_01a_resid <- simulateResiduals(model_fvfm_01a)
plot(model_fvfm_01a_resid)
model_fvfm_01b_resid <- simulateResiduals(model_fvfm_01b)
plot(model_fvfm_01b_resid)
model_fvfm_01c_resid <- simulateResiduals(model_fvfm_01c)
plot(model_fvfm_01c_resid)
summary(model_fvfm_01a)## Family: Gamma ( logit )
## Formula: fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +
## (1 | site) + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## -698.9 -668.6 357.4 -714.9 316
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.03093 0.1759
## year (Intercept) 0.02150 0.1466
## Number of obs: 324, groups: site, 5; year, 14
##
## Dispersion estimate for Gamma family (sigma^2): 0.0214
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.328e-02 9.826e-02 -0.848 0.396709
## urchin_density 8.903e-03 2.766e-03 3.219 0.001287 **
## mean_sand 3.131e-03 1.368e-03 2.288 0.022137 *
## mean_daily_umol 1.607e-05 4.741e-06 3.390 0.000699 ***
## kelp_dry 1.922e-05 4.719e-05 0.407 0.683823
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_fvfm_01b)## Family: Gamma ( logit )
## Formula: fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry +
## site + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## -711.2 -669.7 366.6 -733.2 313
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.01944 0.1394
## Number of obs: 324, groups: year, 14
##
## Dispersion estimate for Gamma family (sigma^2): 0.0212
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.289e-02 7.433e-02 -0.712 0.476774
## urchin_density 9.231e-03 2.757e-03 3.348 0.000813 ***
## mean_sand 3.127e-03 1.368e-03 2.285 0.022321 *
## mean_daily_umol 1.666e-05 4.744e-06 3.512 0.000445 ***
## kelp_dry 2.211e-05 4.683e-05 0.472 0.636884
## sitecarp -2.969e-01 6.435e-02 -4.615 3.93e-06 ***
## siteivee -4.220e-02 6.266e-02 -0.673 0.500657
## sitemohk -1.078e-01 7.375e-02 -1.462 0.143674
## sitenapl 2.468e-01 6.493e-02 3.801 0.000144 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_fvfm_01c)## Family: Gamma ( logit )
## Formula: fvfm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## -618.8 -592.3 316.4 -632.8 317
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 0.02717 0.1648
## Number of obs: 324, groups: year, 14
##
## Dispersion estimate for Gamma family (sigma^2): 0.0288
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.947e-02 6.342e-02 0.938 0.3483
## urchin_density 4.859e-03 3.094e-03 1.570 0.1163
## mean_sand 2.123e-03 1.419e-03 1.496 0.1347
## mean_daily_umol 8.661e-06 4.547e-06 1.905 0.0568 .
## kelp_dry -4.695e-05 5.563e-05 -0.844 0.3987
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_fvfm_01a, model_fvfm_01b, model_fvfm_01c)## df AIC
## model_fvfm_01a 8 -698.8904
## model_fvfm_01b 11 -711.2456
## model_fvfm_01c 7 -618.8032v. max height
ggplot(CWM_fixed, aes(x = urchin_density, y = max_height_cm_fixed)) +
geom_point(aes(col = site)) +
geom_smooth(aes(col = site), method = "lm") +
theme(legend.position = "none") +
facet_wrap(~site)## `geom_smooth()` using formula 'y ~ x'
ggplot(CWM_fixed, aes(x = max_height_cm_fixed)) +
geom_histogram(binwidth = 1)
# full
model_max_height_cm_01a <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|site) + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
# site as fixed instead of random effect
model_max_height_cm_01b <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + site + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
# no site
model_max_height_cm_01c <- glmmTMB(max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol + kelp_dry + (1|year),
data = CWM_fixed, family = gaussian(link = "identity"))
# plot residuals
model_max_height_cm_01a_resid <- simulateResiduals(model_max_height_cm_01a)
plot(model_max_height_cm_01a_resid)
model_max_height_cm_01b_resid <- simulateResiduals(model_max_height_cm_01b)
plot(model_max_height_cm_01b_resid)
model_max_height_cm_01c_resid <- simulateResiduals(model_max_height_cm_01c)
plot(model_max_height_cm_01c_resid)
summary(model_max_height_cm_01a)## Family: gaussian ( identity )
## Formula: max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + (1 | site) + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## 2348.2 2378.5 -1166.1 2332.2 316
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 30.43 5.517
## year (Intercept) 12.02 3.467
## Residual 69.83 8.357
## Number of obs: 324, groups: site, 5; year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 69.8
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.045e+01 2.899e+00 7.055 1.73e-12 ***
## urchin_density 2.239e-01 6.910e-02 3.240 0.0012 **
## mean_sand 1.558e-01 3.744e-02 4.162 3.16e-05 ***
## mean_daily_umol 5.403e-04 1.243e-04 4.347 1.38e-05 ***
## kelp_dry 7.651e-04 1.324e-03 0.578 0.5635
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_max_height_cm_01b)## Family: gaussian ( identity )
## Formula: max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + site + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## 2334.4 2376.0 -1156.2 2312.4 313
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 10.87 3.297
## Residual 68.97 8.305
## Number of obs: 324, groups: year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 69
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 16.9548387 1.9815632 8.556 < 2e-16 ***
## urchin_density 0.2277658 0.0687214 3.314 0.000919 ***
## mean_sand 0.1585724 0.0373990 4.240 2.24e-05 ***
## mean_daily_umol 0.0005753 0.0001225 4.695 2.66e-06 ***
## kelp_dry 0.0008155 0.0013166 0.619 0.535630
## sitecarp -2.6336951 1.8341814 -1.436 0.151032
## siteivee 6.5766697 1.8002934 3.653 0.000259 ***
## sitemohk -0.5999541 2.0136300 -0.298 0.765744
## sitenapl 12.5192031 1.7950750 6.974 3.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model_max_height_cm_01c)## Family: gaussian ( identity )
## Formula: max_height_cm_fixed ~ urchin_density + mean_sand + mean_daily_umol +
## kelp_dry + (1 | year)
## Data: CWM_fixed
##
## AIC BIC logLik deviance df.resid
## 2432.7 2459.1 -1209.3 2418.7 317
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## year (Intercept) 19.52 4.419
## Residual 94.88 9.741
## Number of obs: 324, groups: year, 14
##
## Dispersion estimate for gaussian family (sigma^2): 94.9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.571e+01 1.696e+00 15.160 < 2e-16 ***
## urchin_density 1.504e-01 7.831e-02 1.920 0.05485 .
## mean_sand 1.011e-01 3.872e-02 2.610 0.00904 **
## mean_daily_umol -8.313e-06 1.205e-04 -0.069 0.94501
## kelp_dry -2.262e-04 1.502e-03 -0.151 0.88025
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model_max_height_cm_01a, model_max_height_cm_01b, model_max_height_cm_01c)## df AIC
## model_max_height_cm_01a 8 2348.228
## model_max_height_cm_01b 11 2334.373
## model_max_height_cm_01c 7 2432.673b. specific CWM using FD::dbFD()
i. Bullito
Only one species, skip.
ii. Arroyo Quemado
Only one species, skip.
iii. Isla Vista
# get a site specific trait by species matrix
ivee_tbspp <- specif_tbspp("Isla Vista")
# get a site specific transect by species matrix
ivee_commat <- specif_commat("ivee", ivee_tbspp)
# making sure the rownames in the trait matrix == columns in the community matrix
rownames(ivee_tbspp) == colnames(ivee_commat)## [1] TRUE TRUE TRUE TRUE TRUE TRUE# calculate CWM for site
CWM_ivee <- functcomp(x = ivee_tbspp, a = ivee_commat,
CWM.type = c("all")) %>%
rownames_to_column("sample_ID") %>%
rename_with(.fn = ~ paste0(., "_specif", sep = ""), .cols = 2:ncol(.)) %>%
left_join(., CWM_fixed, by = "sample_ID") %>%
mutate(site = "ivee")iv. Mohawk
# get a site specific trait by species matrix
mohk_tbspp <- specif_tbspp("Mohawk")
# get a site specific transect by species matrix
mohk_commat <- specif_commat("mohk", mohk_tbspp)
# making sure the rownames in the trait matrix == columns in the community matrix
rownames(mohk_tbspp) == colnames(mohk_commat)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# calculate CWM for site
CWM_mohk <- functcomp(x = mohk_tbspp, a = mohk_commat,
CWM.type = c("all")) %>%
rownames_to_column("sample_ID") %>%
rename_with(.fn = ~ paste0(., "_specif", sep = ""), .cols = 2:ncol(.)) %>%
left_join(., CWM_fixed, by = "sample_ID") %>%
mutate(site = "mohk")v. Carpinteria
Only one species, skip.
c. putting all the data frames together
df <- CWM_mohk %>%
select(sample_ID, thickness_mm_specif, tdmc_specif, site) %>%
rbind(., (CWM_ivee %>% select(sample_ID, thickness_mm_specif, tdmc_specif, site))) %>%
rbind(., (CWM_bull %>% select(sample_ID, thickness_mm_specif, tdmc_specif, site))) %>%
left_join(., CWM_fixed, by = "sample_ID")d. variance decomp using cati::traitflex.anova()
tdmc_tf <- traitflex.anova(~ site.x, specif.avg = tdmc_specif, const.avg = tdmc_fixed, data = df)
tdmc_tf##
## Decomposing trait sum of squares into composition turnover
## effect, intraspecific trait variability, and their covariation:
## Turnover Intraspec. Covariation Total
## site.x 0.091056 0.1393 0.11144 0.3418
## Residuals 0.772176 1.4119 -0.52404 1.6600
## Total 0.863232 1.5512 -0.41260 2.0018
##
## Relative contributions:
## Turnover Intraspec. Covariation Total
## site.x 0.04549 0.06959 0.05567 0.1707
## Residuals 0.38574 0.70530 -0.26178 0.8293
## Total 0.43122 0.77489 -0.20611 1.0000
##
## Significance of testable effects:
## Turnover Intraspec. Total
## site.x 0.00016746 0.00064961 4.141e-07plot(tdmc_tf)
thickness_tf <- traitflex.anova(~ site.x, specif.avg = thickness_mm_specif, const.avg = thickness_mm_fixed, data = df)
thickness_tf##
## Decomposing trait sum of squares into composition turnover
## effect, intraspecific trait variability, and their covariation:
## Turnover Intraspec. Covariation Total
## site.x 0.12055 0.26870 0.16142 0.55067
## Residuals 0.56746 0.55187 -0.63954 0.47979
## Total 0.68801 0.82056 -0.47812 1.03046
##
## Relative contributions:
## Turnover Intraspec. Covariation Total
## site.x 0.1170 0.2608 0.1566 0.5344
## Residuals 0.5507 0.5356 -0.6206 0.4656
## Total 0.6677 0.7963 -0.4640 1.0000
##
## Significance of testable effects:
## Turnover Intraspec. Total
## site.x 2.9824e-07 3.6505e-14 8.7061e-27plot(thickness_tf)